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Abstract 

Atmospheric turbulence models are necessary for the design 
of both inlet/engine and flight controls, as well as for studying 
coupling between the propulsion and the vehicle structural 
dynamics for supersonic vehicles. Models based on the 
Kolmogorov spectrum have been previously utilized to model 
atmospheric turbulence. In this paper, a more accurate model is 
developed in its representative fractional order form, typical of 
atmospheric disturbances. This is accomplished by first scaling 
the Kolmogorov spectral to convert them into finite energy von 
Karman forms and then by deriving an explicit fractional 
circuit-filter type analog for this model. This circuit model is 
utilized to develop a generalized formulation in frequency 
domain to approximate the fractional order with the products of 
first order transfer functions, which enables accurate time 
domain simulations. The objective of this work is as follows. 
Given the parameters describing the conditions of atmospheric 
disturbances, and utilizing the derived formulations, directly 
compute the transfer function poles and zeros describing these 
disturbances for acoustic velocity, temperature, pressure, and 
density. Time domain simulations of representative 
atmospheric turbulence can then be developed by utilizing these 
computed transfer functions together with the disturbance 
frequencies of interest. 

Introduction 

This paper addresses the need for a model that simulates 
atmospheric disturbance, in both the time and frequency 
domain, over a wide range of altitudes and variations in 
atmospheric turbulence conditions, that is relatively easy to 
implement and representative of the actual fractional order 
nature of atmospheric turbulence. This is applicable to both 
propulsion system flow field type disturbances as well as 
vehicle gust loads. 

Atmospheric turbulence have been studied for some time, 
especially in the field of Atmospheric Sciences, Nastrom 
(1985), and Fairall (1991), and various models have been 
developed. These models are primarily based on the so called 
Kolmogorov spectrum, originally developed by Tatarski (1961), 
based partly on the studies of turbulence by the Russian 
mathematician Andrei Kolmogorov (Kolmogorov (1941) (a) and 
(b)). The Kolmogorov spectrum has an energy level that 
approaches infinity as frequency approaches zero. This 


characteristic makes it difficult to implement these types of 
models in the time domain. A suitable approximation to the 
Kolmogorov model that is commonly used with a finite energy 
spectrum is the von Karman type model, originally referred to 
as the isotropic-turbulence spectrum, (Houbolt (1964)). 
However, simulating the von Karman type models in the time 
domain is still problematic because of their fractional order. 
The Kolmogorov model has also been extended (Tank (1994)) 
to develop a baseline of atmospheric turbulence for the High 
Speed Civil Transport (HSCT). The Tank model also covers 
atmospheric acoustic wave disturbance modeling utilizing the 
von Karman spectral. Hoblit (1988) introduced the Dry den 
model approximation to the fractional order von Karman 
atmospheric model. But this model is second order compared to 
the 5/3 fractional order of the acoustic velocity atmospheric 
turbulence spectral. Thus, the Dryden model underestimates the 
atmospheric disturbances, increasingly with frequency. 

To alleviate some of these difficulties, the Tank model for the 
von Karman approximation is utilized in this paper to derive an 
explicit electrical circuit analog of atmospheric disturbances. 
This circuit analog will act as a low pass filter. The circuit 
elements are explicitly computed as functions of atmospheric 
parameters, such as eddy dissipation rate and integral length 
scale. However, like the actual atmospheric disturbances, 
Nastrom (1985), the circuit order also turns out to be fractional 
order, which makes it difficult to simulate. Thus, the circuit 
model is used as the basis in this development to derive integer 
order transfer function (TF) approximations to the fractional 
order model. These TF approximations are a product of first order 
poles and zeros, which are determined as a function of the 
parameters describing an atmospheric disturbance. This approach 
alleviates the manual process of hand fitting the approximation 
every time an atmospheric parameter is changed. 

Fairall (1991), Tank ((1994) and (1996)) and others approach 
atmospheric disturbance modeling probabilistically, as an 
exceedance for controls design purposes. That means that the 
probability of a time to failure is computed for controls design. This 
time to failure is associated with a controls design that can tolerate 
a maximum specified disturbance, and computes the probability, 
in terms of flight miles or hours, that an atmospheric disturbance 
will exceed this threshold. In this paper the approach will be to 
compute the worst case disturbance expected, and assume that the 
control system will be designed to handle this disturbance. Thus, 
this paper will not cover exceedance rates , which can be computed 
separately, based on the worst disturbance expected and the 
specifics of the controls design. 
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Kopasakis (2010) extends the work discussed in this paper to 
cover considerations for atmospheric turbulence specifications 
for supersonic propulsion systems. This work also provides an 
example for a supersonic inlet shock position controls design in 
the presence of atmospheric turbulence. 

This paper is organized as follows. First, the Kolmogorov 
form of the atmospheric disturbance spectrum is presented 
followed by the Tank model for the von Karman spectrum 
forms of the acoustic disturbances. This is followed with 
formulations of the equivalent fractional order TF 
approximations of atmospheric turbulence. Finally, time 
domain atmospheric disturbances are discussed for the derived 
TF approximations, followed by concluding remarks. The 
formulations covered in the body of the paper are based on an 
electrical circuit analog developed for different types of 
atmospheric disturbances, which is covered in Appendix B. 
While the detailed derivations of these formulations are covered 
in Appendix C. 


Kolmogorov Form of the Atmospheric 
Disturbance Model 

Tank ((1994) and (1996)) utilizes a Kolmogorov one- 
dimensional locally isotropic atmospheric turbulence spectrum, 
mathematically developed in Tatarski (1961), which represents 
the spectral density of a structured random field of atmospheric 
turbulence as 


S t (k) = a t £ 2/3 k~ 513 


(i) 


exceeded. The quantity v is the kinematic viscosity, and v\ is a 
velocity fluctuation in an atmospheric region of size /. 

At lower altitudes, the eddy dissipation rate can vary 
significantly with altitude and severity of turbulence (Fairall 
(1991)) as well as globally (Nastrom (1985)). Tank (1994) used 
a worst value of 8.6e-5 based on data collected at altitudes 
ranging from 25,000 to 40,000 ft (~ 7.5 to 12 km), which is 
about four times the appropriate value of 8 for North Atlantic 
cruise altitudes. For supersonic cruise altitude, around 60,000 ft 
(-18 km), the same value of e can be used, as s is relatively 
constant at altitudes above 20,000 ft (-6 km). 

According to Tatarski (1961), Equation (1) holds for 
turbulence lengths much greater than what would be considered 
a microstructure (10’s of meters) for eddy dissipation and much 
less than the long length scale (that is considered to be about 
400 km). Some sparse measurements performed by Obukhov 
(1949) in the lower troposphere agree with Equation (1), as well 
as other measurements by Nastrom (1985). The constant a t is 
constant for each type of disturbance, given by Tank 1994, to 
fit observed data: 

a /= 0.15 

a v = 0.2 

qlt= 0.39 

clp= 0m05(Po/T o ) 2 


(longitudinal wind velocity gust, 
dimensionless) 

(vertical or horizontal wind velocity 
gust, dimensionless) 

(temperature disturbance, K 2 s 2 m 2 ) 
(pressure disturbance, hPaVm 2 - 
Tank 1994, multiply a p by 100 to 
change units to Pa 2 s 2 m 2 ) 


In this equation the quantity 8 signifies the eddy dissipation rate, 
in units of (energy/(mass x time)) = (m 2 /sec 3 ), and k signifies 
the wavenumber, in units of (rad/m) or (cycles/m). The 
subscript t in the atmospheric turbulence spectral density, S t (k), 
signifies the type of disturbance. Based on these definitions, the 
units of S t (k) are (m/sec) 2 /(rad/m) (Tatarski (1961)), where the 
units of rad or Hz are not affected by raising them to a power. 
Or S t (k ) will have units of (m/sec) 2 /(cycle/m) (Tank (1994)), as 
rad = cycle/27i, with the 271 factor that multiplies the spectral 
density for this conversion absorbed into the constants a t terms. 
In this treatment, more convenient units in terms of Hz will be 
utilized as ((m/sec) 3 /Hz), by substituting cycles with Hz* sec, 
which allows displaying results in terms of the acoustic wave 
velocity in (m/sec) /Hz. 

Tatarski (1961) includes a more detailed treatment of 8 along 
with a detailed derivation of Equation (1). In brief, 

8 = v * V 2 // 2 

1 ' is called the eddy dissipation rate and it 
represents the energy dissipated as heat per unit mass per unit 
time due to a velocity fluctuation that occurs in an atmospheric 
region of size (length) /, due to a flow instability that takes place 
when a certain critical Reynolds number value is exceeded 
(Tatarski (1961)). This critical Reynolds number cannot be 
determined precisely, but for worst case atmospheric turbulence 
calculations it would be correct to assume that this number is 


Note Tank 1994 does not define the symbols P Q and T 0 for 
pressure and temperature in the expression for a p. McMinn 
1997, who is referencing Tank 1994, defines these symbols as 
standard atmospheric quantities. However, the assumption here 
is that Tank 1994 meant to define these quantities as freestream 
atmospheric conditions. The reasons for this assumption are as 
follows: a) because if these quantities were referring to standard 
atmospheric conditions, a p would turn out to be a constant, just 
like the rest of the a t constants (no need to display it in a variable 
form), and b) because Tank 1994 describes the atmospheric 
turbulence measurements taken for pressure as being freestream 
measurements, and as such, it would have not made much sense 
to convert these measurements into stagnation quantities, either. 

For a flight vehicle encountering an atmospheric disturbance, 
the disturbance frequency is defined as 


r Ma 1 a /T 

f = = kMa 

X 


( 2 ) 


Where M is the vehicle Mach number, a is the local speed of 
sound (m/sec), X is the wave length of the atmospheric 
disturbance (m/cycle), and k is 1/A, (cycles/m). The disturbance 
frequency is based on an aero vehicle traveling at a certain 
altitude and Mach number that encounters an atmospheric 
disturbance with a certain wavelength. There would be a 
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corresponding frequency based on Equation (2) around which 
atmospheric disturbances will be exhibited. However, a control 
design would need to be able to sufficiently attenuate these 
disturbances under a wide range of flight and atmospheric 
conditions. Thus, accounting for the full operating envelope, a 
controls design would need to attenuate atmospheric 
disturbances in a wide range of frequencies. Normally, lower 
frequency atmospheric disturbances will be easier for the 
control design to handle. 

The speed of sound is related with temperature as 

a = jy RT s ( 3 ) 


(m 3 /sec 2 )/rad and temperature is in (K 2 m/rad). As shown in 
this figure, at the long length scale (i.e., beyond 400 km) 
the power of turbulence decreases as the -3 power with 
decreasing wavelength (increasing frequency), thereafter 
decreasing as the -5/3 power which agrees with 
Equation (1). Therefore, Equation (1) shows a good 
representation of the atmospheric spectral density with a 
wavelength smaller than 400 km. 

The Tank Model of the Atmospheric 
Disturbance 


where T s is the static temperature in K. 

According to the International Standard Atmosphere (ISA) or 
(Anderson (2000)), the temperature decreases at a rate of 
approximately 6.5 K per km, from sea level up to 36,000 ft 
(-11 km). From 36,000 to 65,000 ft (-11 to 20 km) which 
would be the approximate cmise altitude of a supersonic 
transport, the temperature remains constant at about 216 K. 
Thus, temperature (K) as a function of altitude, h (km), can be 
expressed as 


As an alternative to the Kolmogorov spectmm, Tank (Soreide 
and Tank (1996) and (1997)) scaled the von Karman spectmm 
to fit the Kolmogorov model in the limit (i.e., for large k ), which 
compares well with other data. For longitudinal disturbances 
this spectmm is 


S UVK {k)=2.1& 1 ^L^^ 


2 

1 + (1.339(2tt )Lk) 2 5/6 


( 5 ) 


T s (h) = T oh - 6 . 5 ( 7 ? - h oh ) 


(4) For transverse disturbances this spectmm is 


where T 0 h is the temperature at an arbitrary reference height h 0 h, 
which is accurate up to 36,000 ft and then stays constant above 
that. Based on that, at supersonic cmise conditions the speed of 
sound would be approximately 295 m/sec. 

Figure 1 shows the zonal (east- west) winds, meridional (north- 
south) wind, and potential temperature power spectral as reported 
in Nastrom (1985). The units of the wind power spectral are in 



Figure 1. — Wind and Potential Temperature Spectra as 
reported by Nastrom (1985). Note: for clarity, the 
meridional wind and potential temperature spectra have 
been shifted one and two decades to the right, 
respectively. 


S v , VK {k)=2Js 2/3 L 5/ 3 


l + -(l .339(2 ji )Lkf 


l + (l.339(2n)Lkf 


11/6 


( 6 ) 


The main difference of this von Karman form (compared to 
Kolmogorov model) is that the disturbance spectral levels off at 
low frequencies. Another difference is that in the von Karman 
form of the Tank model, the integral length scale, Z, is explicitly 



Frequency, Hz 

Figure 2. — Acoustic wave velocity spectral comparisons for the 
Kolmogorov and von Karman spectral. 
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employed, Equations (5) to ( 6 ). The integral length scale is 
related to the outer scale length, L 0 (the length of the 
atmospheric turbulence patch) by L=L 0 1 14.7 for the 
longitudinal disturbance, and L=2L 0 /19 for the transverse 
disturbance (Soreide and Tank (1996)). Figure 2 shows a plot 
of the Kolmogorov acoustic wave velocity spectral in 
(m/sec/Hz) compared to the corresponding von Karman 
spectral for the longitudinal and transverse acoustic wave 
velocity disturbances. Different values of L will produce the 
same curves with either a higher or lower, low-frequency 
asymptote as shown in Figure 3 for the longitudinal 
disturbance. For comparison, a value of L=1 62 m pertains to an 



Frequency, Hz 

Figure 3. — Longitudinal acoustic wave velocity spectral 
comparisons between the Kolmogorov and the von Karman 
forms of four different integral scale lengths. 



Figure 4. — Longitudinal acoustic wave velocity spectral 
comparisons for different integral scale lengths. 


atmospheric turbulence patch of approximately 1 1 km for the 
longitudinal disturbance, while a value of L =10 km pertains to 
a turbulence patch of approximately 147 km long. Tank used a 
value of L=1 62 m, which he called as standard in the airplane 
industry. Since a typical control design can sufficiently 
attenuate disturbance frequencies well beyond 1 Hz, based on 
Figure 3 differences in the value of L and thereby, differences 
in the lower frequency asymptote, will have negligible effect on 
the control design. For completeness, the Kolmogorov 
spectrum for the longitudinal acoustic wave velocity is shown 
in Figure 4 for different values of eddy dissipation rates. As can 
be deduced from Figure 4, orders of magnitude difference in 
eddy dissipation rate does not produce a corresponding 
appreciable change in the Kolmogorov spectral density. For 
instance, approximately four orders difference in s only makes 
approximately a factor of 8 difference in velocity as 


W) = 


/ \ 2/9 

c 2 


V £ 1 ) 


Sitf) , ((m/sec)/Hz) 


( 7 ) 


Where Si is a known spectral density with an eddy dissipation 
rate Si, and & is the calculated spectra density for a different 
eddy dissipation rate 82 . 


Fractional Order Fit of Atmospheric 
Turbulence Model 

Atmospheric turbulence, as shown in Figure 1, and as 
described by Kolmogorov and the Tank models in Equation (1) 
and Equations (5) and ( 6 ), is fractional order. In Appendix B a 
circuit analog of the Tank model von Karman form is utilized, 
which serves as the basis for deriving integer order pole-zero 
product TF approximations to the fractional order atmospheric 
disturbances. The reason for the need to derive approximations 
to the fractional order equations is because of the difficulty of 
explicitly or numerically solving fractional order differential 
equations. The reason is fractional order derivatives, unlike 
integer order derivatives, do not obey the locality law (i.e., the 
limit theorem), as further explained in Appendix B. 

The idea behind the integer order TF approximation for a 
fractional order TF is as follows (see Fig. C.l): Starting at a 
frequency near the beginning of the equivalent 3 dB point of the 
fractional order TF, an integer order TF approximation can be 
developed symmetrically centered about the fractional order 
TF, like a descending staircase shape, by interleafmg integer 
order poles and zeros. As the number of steps of this staircase 
TF approximation increase, the steps become shorter, 
eventually collapsing to the straight line of the fractional order 
TF as the number of poles and zeros of this approximation is 
increased to infinity. For this approach to work, the frequency 
of the poles and zeros need to be related to the atmospheric 
disturbance parameters and also be derived in a way such that 
the staircase TF approximation is symmetrically centered about 
the fractional order TF. 
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Based on the fractional order circuit analog of Equations (5) 
and (6), see Appendix B for detailed derivations, the values for 
the equivalent capacitance and resistance of atmospheric 
turbulences are determined as 


Q=- 


(^8 2/3 ) 1/x (27iMa) 

R t =13390(2n)(a t £ 2/3 f x L 


and the corresponding natural frequency is computed as 


00 „ 


RfC t 


( 8 ) 


(9) 


( 10 ) 


where the adjustment factor K m is 1 , but represented in this form 
in case any adjustments need to be made to the atmospheric 
disturbance natural frequency. Note that the natural frequency, 
co n , does not depend on the different type of disturbances or the 
eddy dissipation rate, as these factors cancel out when 
resistance and capacitance are multiplied together in 
Equation (10). 

Utilizing these circuit analog parameters, a time domain 
disturbance can be derived (see Appendix C for detailed 
derivations) based on an integer order TF approximation for 
atmospheric turbulence, given the atmospheric disturbance 
parameters 8, Z, q , t , and r for the units conversion factor (with 
q=xr ; x=5/3 and 7=1/3 for acoustic disturbances and l A for 
temperature and pressure) as 


m z 

ru/®zi +0) 

W t ,o = w, (1 1) 

1 1 (v®/"’ + 1) 

1 

where W t is series of input sinusoids with unit amplitude 
frequency components distributed over the frequency range of 
interest, which is discussed more in Kopasakis (2010) for 
propulsion system type disturbances. The frequencies of the 
poles can be computed as 


K<x>pi^Hpi J[ J[ i^Hpi /® pi-j + l) 

(dpi = i = 2,3,...,m p (12) 

1 0 nf2/-i)£/ jq (c,) [[pi /o) zi _j + 1)- 1 
j = 1 

where the first pole is computed as follows 

(dpi =co„(l0 1 i9-lU (13) 

and the frequencies of the zeros can be computed as 


Kazi^HziYYi^Hzi/^zi-j + 1 ) 

o »zi= / — i' = l,2,...,m z (14) 

lO-^n^/ffl^ + lJ-l 

7=1 

For n number of frequency decades desired to be estimated, the 
number of poles and zeros in the TF approximation can be 
determined as 


m p =(n- l)/r| (15) 

m z = m p - 1 (16) 


For a desirable pole-zero density pair, p pz , per decade to be used 
to approximate the fractional order disturbance 


1 

pz 


(17) 


and the terms co Hp i and oo Hz i in Equations (11) and (12) can be 
computed as 


(d Hpi = (d n (l0d^) -if \ i = 2,...,m p (18) 


(dHzi = (d„ (l 0 2t19 ' - lj /? , i = l,2,...,l« z (19) 

The utility of co H P i and co Hzi are to maintain symmetry, so that 
the staircase pole-zero approximation is symmetrically located 
on top of the fractional order TF asymptote. The proportional 
gains Ka pi and K mi in Equations (12) and (14) are reserved for 
final adjustments that may be needed to these frequencies for 
more closely approximating the fractional order TFs 
representing atmospheric disturbance. 


Longitudinal and Transverse Acoustic Wave 
Turbulence 

The longitudinal and transverse acoustic wave atmospheric 
disturbances are in the form of pure wind gusts in the axial 
direction of vehicle motion and in the vertical direction of 
motion respectively. Based on the fractional order fits 
determined in Appendix C, using Equations (8) to (19), with 
q=xr= 5/9 (r=l/3), then K t> m in Equation (11), based on 
inspection of Equations (5) and (6), is 

X />fit =(5.4e 2/3 Z? /3 )f /3 (20) 

for the longitudinal disturbance, and 

^v,fi t =(2.7s 2/3 ^3) 1/3 (21) 
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for the transverse disturbance. For n = 3 (i.e., the TF fit over the 
span of 3 decades in frequency), the proportionality factor 
adjustments to improve these fits have been found to be 
(Appendix C), 

*/,« — [^co?z ? , K^ z f ] (22) 

= [2.4; 1 1 1/2.4 1/1.5; 1 1 l] 


a good job in approximating these disturbances. More accuracy 
can be achieved around the knee of these spectral by increasing 
the density of pole-zero pairs per frequency decade in this 
region. A certain inaccuracy at the lower frequencies would be 
acceptable for a feedback control design, since a typical control 
design should have no problem attenuating disturbances at this 
low frequency range. 


for the longitudinal acoustic wave, and 

^V,G> = [^0072 ? Kffipi , K c oz i ] (23) 

= [4.27; 1 1 1/2.4 1/1.5; 1 1 l] 

for the transverse. These TF fits are compared against the Tank 
von Karman spectral of Equations (5) and (6), with the units 
conversion exponent r, as shown in Figures 5 and 6 for s=8.6e- 
5 (m 2 /sec 3 ) andZ=762 m. As shown in these figures, the fits do 



Figure 5. — Longitudinal von Karman spectral, final adjusted 
TF fit (e = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Figure 6. — Transverse von Karman spectral, final adjusted 
TF fit (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


Potential Temperature Turbulence 

For the sake of simplicity, the terminology for temperature 
turbulence/disturbances, T, is utilized in this development to 
indicate potential temperature. Maintaining potential temperature 
in the formulations and plots versus actual temperature, also 
allows for direct comparisons between Kolmogorov (see Fig. 1), 
von Karman, and fractional order fit models. Potential 
temperature is defined as the temperature that a parcel of fluid at 
pressure, P , would acquire if adiabatically brought to a standard 
reference pressure P sh usually 100,000 Pa. To convert to actual 
static temperature, T s , the levels displayed in the equations and 
the graphs should be adjusted by the following equation: 

Ts = T (I_P’ (24) 

where R is the universal gas constant and c p is the specific heat 
capacity at constant pressure. R/c p = 0.286 for air. Notice, that at 
altitudes near sea level, potential and actual temperatures are 
about the same, while they increasingly deviate with altitude (eg. 
T, -0.57 at 50,000 feet). 

Atmospheric temperature turbulence causes both temperature 
as well as acoustic velocity disturbances due to the change in the 
speed of sound. For the most part, temperature disturbances result 
in vertical displacement of air parcels (so called gravity waves). 
Therefore, acoustic velocity disturbances caused by temperature 
will generate vertical wind gusts that add with any transverse 
acoustic velocity gusts. However, for the propulsion system, 
according to Ashun (2004), the wing forebody will turn a vertical 
gust into a longitudinal gust, by multiplying the vertical gust by 

the conversion factor ( M ^ - 1)/ ^/m 2 -1 . This factor amounts 
to 0.63 for Moo = 2.35, and can be ignored for worst case 
purposes. However, this conversion factor decreases with speed 
and can be taken into account, especially at lower supersonic 
speeds, for longitudinal acoustic gusts due to temperature 
turbulence. In Appendix B, Section B.2, first the Kolmogorov 
spectrum of the longitudinal acoustic wave and that of the 
temperature, based on Equation (1) are plotted (Fig. B.6). This 
results in parallel spectral lines with frequency, similar to the 
Kolmogorov spectra shown in Figure 2. Then the von Karman 
temperature spectral is constructed by scaling the horizontal von 
Karman type acoustic wave spectral (i.e., scaling it by the 
difference in magnitudes of these Kolmogorov spectra as shown 
in Fig. B.6). 

Temperature turbulence spectral density, as can be seen in 
Figure 1, also follows the 5/3 law, but its units are in terms of 


NAS A/TM— 20 1 0-2 1 696 1 /REV 1 


6 


Kelvin squared. Thus, in order to convert the units to Kelvin, 
the exponent r becomes V2, which makes the fractional order 
q= 5/6. Based on this, the TF fit performed for temperature is 
done the same way as with acoustic wave in the previous 
section (i.e., utilizing Equations (8) to (19), with the additional 
relations. 


This relation can be applied to Equation (25) to compute 
Kr,fit(acoustic) for the acoustic wave velocity disturbance due to an 
atmospheric temperature fluctuation as 

^r,fit(acousic) = ^140^3^3 (28) 


=Vl4.0s 2/3 Z 5/3 (25) 

Kt,B) — hoon > K a pi , K az i ] ( 26 ) 

= [1.5; 1 1 1/1.1 1/1.2; 1 1 1] 


Figure 7 shows this TF fit and the spectral of Equation (5) (i.e., 
Equation (5) substituting the scaled Equation (25) for the 
numerator and by also applying the units conversion factor in 
the denominator, for s=8.6e-5 an L=1 62 m. Note that the units 
of temperature include (m/sec) 1/2 , which is preserved in order to 
maintain consistency with the derivations. However, the actual 
units would not include the V 2 power factor. Also, as discussed 
previously, temperature in this development refers to potential 
temperature. 

The proportionality factors, K t , for all the fits can also be 
considered part of the formulations. Their accuracy at higher 
frequencies starts deviating somewhat with higher values of L 
(integral length scales of few thousand meters and above), 
values that are considered to be outside the typical range. As 
discussed before, L=162 m is considered standard in airplane 
industry. 

Temperature turbulence also causes acoustic wave gusts due 
to change of the local speed of sound. By perturbing the relation 
of the speed of sound and temperature in Equation (3), and 
substituting the resulting expression for A a into the Mach 
number equation of M=v/a, the relation between the change in 
temperature and acoustic velocity can be obtained as 


Av = 


MyR 

2a 0 


AT 


(27) 



10-3 io-2 10-I 10O 101 102 

Frequency, Hz 

Figure 7. — Temperature von Karman spectral and 
its TF fit (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


Equation (28) is applicable to transverse or vertical wind gust 
disturbances generated by temperature turbulence, and accurate 
in the worst case sense for longitudinal wind gusts, especially 
at higher supersonic Mach numbers. Generally, longitudinal 
wind gusts due to temperature variations can be expressed as 

follows by utilizing the conversion factor (M^ - 1)/ -1 

Kt ,fit(longitudinal) = (29 ) 

2a 0 yjM^ -1 

This TF fit for the acoustic velocity disturbance at supersonic 
cruise altitude, at Mach 2.35, for L=762 m and two different 
values of c is shown in Figure 8. The TF fit for c=8.6xl0“ 5 
(m 2 /sec 3 ) and for two different values of integral scale lengths is 
shown in Figure 9. As can be seen from Figure 8 or Figure 9 
compared to Figure 5, the acoustic velocity gusts produced by 
fluctuations in temperature are much higher in amplitude than 
those produced by pure acoustic velocity gusts. 

As shown in Appendix B, Figure B.8, for the combined 
acoustic wave, the acoustic wave velocity due to temperature 
gust dominates in the low frequency range and at higher 
frequencies both acoustic waves affect the total. At worst case, 
it can be assumed that both effects of the acoustic wave velocity 
can combine. It is also assumed that both the longitudinal and 
transverse acoustic wave velocities will not combine to produce 
worst case affects. For engines under the wing, the transverse 
wave is assumed to be converted to a longitudinal wave via the 
wing forebody (Ahsun (2004)), see Appendix B. 



Figure 8. — Von Karman acoustic wave velocity spectral 
due to temperature gust and its TF fits for different 
eddy dissipation rates (L = 762 m). 
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Figure 9. — Von Karman acoustic wave velocity spectral due to 
temperature gust and its TF fits for different integral scale 
lengths (e = 8.6e-5 m 2 /sec 3 ). 


Pressure Turbulence 

The longitudinal velocity wave spectral and the pressure 
spectral of Equation (1) was utilized the same way as described 
in the previous section for the temperature to scale the von 
Karman type acoustic wave spectral to come up with the 
pressure spectral of atmospheric disturbances, see Section B.3. 
The units conversion exponent for pressure is the same as that 
for temperature (i.e., r=l/2). Utilizing Equation (B.25), with 
this unit conversion factor (same fractional order as 
temperature, q= 5/6), K p> f lt can be expressed as 

K P , fit = (30) 


and K Pi o 0 was computed as follows: 

~ 0 pi , K^ z i ] ( 31 ) 

= [l.5; 1 1 1/1.1 1/1.2; 1 1 l] 

Figure 10 shows this TF fit for c=8.6e-5 (m 2 /sec 3 ) andZ=762 m 
compared to the scaled von Karman type spectral — i.e., 
Equation (30) substituted for the numerator and proportional 
factor of Equation (5) and the denominator of Equation (5), 
raised to the power r. 

Density Turbulence 



Figure 10. — Pressure von Karman spectral and its TF fit 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Figure 1 1 . — Density von Karman spectral and its TF fit 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


density disturbance with temperature and pressure can be 
derived as (by perturbing the state equation) 


P +AP P 

A = £*_ 

R(T 0 +M) RT 0 


(32) 


A plot of density disturbance utilizing Equation (32) is shown 
in Figure 11. 


Simplified Atmospheric Turbulence Model for 
Controls Design 


Density disturbances are contained within temperature and 
pressure fluctuations and its inclusion in a simulation, along 
with pressure and temperature will produce unnecessary 
additional disturbances. However, if need be, the relation of 


The attempt in this section will be to further simplify the 
fractional order TF fit developed in the previous sections, by 
identifying certain characteristics of this model that are 
important for controls design to handle atmospheric turbulence. 
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Based on observations, utilizing Figure 8, differences in eddy 
dissipation rates, 8, produce spectral densities with the same 
frequencies, offset by a fixed magnitude. This offset is 
accounted by differences in TF gain. Calculating these 
frequencies using the formulations developed so far also 
supports this observation. Differences in the fractional order of 
the disturbances, such as longitudinal acoustic wave velocity 
vs. temperature, obviously, will produce different spectral 
frequencies. Unlike differences in eddy dissipation rates that 
produce the same frequencies with different TF magnitudes, 
inspection of Figure 9 shows that spectral with different integral 
length scales, L , will be represented by different frequencies. 
However, different integral length scales only affect the low 
frequency spectral, where a typical feedback controls design 
should have no difficulty attenuating disturbances at such a low 
frequency range. 

Figure 12 shows spectral density TF fits, first, with a fixed 
integral length scale and two different values of 8. This shows 
that the differences in 8 produce a fixed offset in magnitude, but 
with the same pole-zero frequencies. Then Figure 12 also shows 
comparison between two spectral densities with the same value 
of 8 and two different values of L. This later comparison, with 
different values of L, is represented with different pole-zero 
frequencies. However, their high frequency asymptotes match, 
which means that for the purpose of controls design, either of 
these TF fits can be utilized, despite their differences in pole- 
zero frequencies. 



Figure 12. — von Karman longitudinal acoustic wave velocity 
spectral and its TF Fits due to different values of eddy 
dissipation rates and Integral Length Scales. 


Therefore, the frequency domain atmospheric turbulence 
model can be simplified by utilizing the model developed in the 
previous section to calculate the fixed frequencies due to a fixed 
integral length scale and variable eddy dissipation rates. Based 
on that, the simplified model developed for atmospheric 
turbulence, which ignores differences in the pole-zero 
frequencies due to different integral length scales can be 
represented as 


Gla(s) = 70s 2/9 


(.s7 9.2 + 1)(, s7 5 5.0 + 1)(,s7 335.5 + 1) 

(,v / 1 .46 + 1 )(,s730. 1 + l)(,v / 85.7 + l)(,v / 1 593. 1 + 1) 


G^(,) = 56^ Qs7 9.2 + l)(s/ 55.0 + I)(.s7 335.5 + 1) 

(^/l .46 + 1)(7 / 30. 1 + 1)(7 / 85.7 + l)(s / 1 593. 1 + 1) 

q ( .) 943s 2/6 (7/ 33.0 + !)(.?/ 45.6 + !)(■?/ 602,4 + 1) 

( 5 /I.I + l)(s / 25. 1 + 1)(7 / 109.8 + l)(s /8 16.3 + 1) 

q {s) 1?2c? ,/ fi (37-1) MyR (5/33.0 + l)(5/ 45.6 + l)(5/ 602,4 + 1) 

TLAK ’ 2_i a 0 ( 5 /I.I + 1)(5 / 25 . 1 + 1)(5 / 1 09 .8 + l)(s / 8 1 6.3 + 1) 

G ( 5 ) 172 C ?/fi MY * (5/33.0 + l)(5/45.6 + l)(5/602.4 + l) 

TVA ’ a 0 ( 5 /I.I + 1)(5 / 25 . 1 + 1)(5 / 109.8 + 1)(5 / 8 16.3 + 1) 


G P (s ) = 586 (a P / a^s 


2/ 6 (5/33.0 + l)(5/45.6 + l)(5/602.4 + l) 

(5/1.1 + 1)(5/25.1 + 1)(5/109.8 + 1)(5/816.3 + 1) 


(33) 

(34) 

(35) 

(36) 

(37) 

(38) 


Where Gla, Gva, Gt, Gtla , Gtva and Gp are the simplified 
atmospheric disturbance TF for longitudinal and transverse 
acoustic wave velocity, temperature, longitudinal and vertical 
acoustic velocities due to temperature, and pressure 
respectively. 


The temperature related TFs of Equations (35) to (37) are 
expressed in terms of potential temperature, in order to allow 
for direct comparisons for the simplified TF fits. However, in 
order to convert them to actual temperature and associated wind 
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velocities, these TFs should be multiplied by the corresponding 

/ p \ R / C p 

conversion factor of Equation (24) or I — ) 

' P St ' 

The factor (m^ -1 in Equation (36) can be set to 

one for worst case turbulence conditions, especially for 
relatively high Mach numbers, in which case it makes Equation 
(36) and Equation (37) the same. 

The TF fit described in the previous section was carried up to 
200 Hz, which should make it accurate even beyond this 
frequency range. However, the model simplification 
represented by Equations (33) to (38) introduces additional 
error at high frequencies and its validity should be limited to 
about 200 Hz. The frequencies of this model are based on an 
integral length scale of L=7 62 m. This means that this 
simplified model is not unique, and other equivalent models can 
be calculated using different values of L with the equations 
developed in the previous sections. Thus, for a controls design, 
which can adequately attenuate 1 Hz disturbances, this model 
is accurate for values of L as low as 100 m, which is well below 
the expected range of L. A value of L=162 m is considered 
standard in the airplane industry as mentioned before. 

The frequencies calculated in this simplified model are also 
dependent on the vehicle Mach number, see Equations (8) and 
(10). For this simplified model, the Mach number used is 2.35, 
which makes this model sufficiently accurate down to a Mach 
number of about 1.5. However, a simplified model can be 
calculated at any atmospheric Mach number, using the 
equations developed in this paper and the assumptions utilized. 

Time Domain Disturbances 

The TF fits in the previous section can be thought as filters, 
which take input disturbances at desirable frequencies and 
convert them to representative free stream atmospheric 
disturbances. The TF fits contain the representative magnitudes 
of these disturbances. Thus, the discrete time domain 
frequencies, W t , only need to have unity RMS values. For 
instance, a time domain input can be constructed as a sum of 
unit amplitude sinusoids, starting with a low frequency, 
somewhere at the low frequency asymptote of the spectral, and 
continuing up with sinusoids that are dispersed along the high 
frequency asymptote up to frequencies that covers the 
controller bandwidth. Also, the control design can be checked 
with single frequency sinusoids, sequentially simulated. Of 
course, the total energy or the maximum amplitude of the 
disturbance will need to be reasonable. For instance, it may be 
considered that total wind gust will not exceed an amplitude of 
180 mph, or some other reasonable number that depends on 
altitude. For more information on propulsion type disturbance 
considerations, see Kopasakis (2010). 

For most of the spectral density simulations covered so far, 
the worst case eddy dissipation rate for low to moderate 
turbulence was used for altitudes of about 36,000 ft (~1 1 km) 
and above. However, the highest atmospheric turbulence ever 


recorded had an eddy dissipation rate of 1.7e-3 (m 2 /sec 3 ) (Tank 
(1994)). Theoretically, it is possible that even more severe 
turbulence can be encountered according to McMinn (1997). 
McMinn provides tables for altitudes, turbulence severity, and 
associated eddy dissipation rates. Kopasakis (2010) covers 
more in this area of simulating atmospheric disturbances, 
including considerations for atmospheric turbulence 
specifications. Fairall (1991) also covers implementation of 
time domain sinusoidal input disturbances by utilizing Fourier 
sinusoids. 

As mentioned before, previous works cover atmospheric 
turbulence probabilistically, as an exceedance for controls 
design purposes. In terms of probabilistic atmospheric 
turbulence models, there are no comparisons drawn here 
between those models and the non-probabilistic approach 
described in this paper, even though exceedance can also be 
incorporated with this approach. The modeling approach 
developed in this paper is geared towards implementing 
controls design by considering worst case atmospheric 
turbulence. This also allows implementation of control designs 
to cover worst case turbulence encounters, without penalizing 
performance under normal atmospheric conditions, by 
incorporating control scheduling in the design. 

Figure 13 shows a time domain longitudinal acoustic wave 
disturbance utilizing the acoustic wave TF fit discussed in the 
previous section. The input to this TF is a combination of unity 
amplitude sinusoids distributed along the frequency range 
shown in Figure 5. For a single unity amplitude sinusoid at low 
frequency, the amplitude of the acoustic disturbance would be 
approximately 9 m/sec according to Figure 5. For this 
summation of unity amplitude input sinusoids, their amplitudes 
combine to give approximately 15 m/sec maximum 
disturbance. The same approach can be used to construct other 
atmospheric disturbances for time domain simulations. 



Figure 13. — Longitudinal acoustic wave velocity disturbance 
created using the TF fit of the longitudinal acoustic wave 
with unity input sinusoids. 
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In Kopasakis (2010), considerations for atmospheric 
turbulence specifications are provided, with more discussion on 
the frequencies of these turbulences, also involving the vehicle 
Mach number. In this reference a control design example is also 
provided for the supersonic inlet shock control problem, in the 
presence of atmospheric turbulence. 

Conclusion 

In this paper derivations were carried out to approximate the 
fractional order atmospheric turbulence with integer order pole- 
zero transfer function fits for more accurate time domain 
simulations. An existing von Karman type model form of the 
Kolmogorov spectral is utilized initially for the acoustic 


velocity gusts disturbances. This model is scaled in this paper 
in order to also develop the von Karman model forms for other 
type of atmospheric disturbances like temperature, pressure, 
and density. Because of the fractional order of these models, a 
circuit equivalent is developed that is used as the basis to derive 
the integer order pole-zero approximations. Utilizing the 
formulations presented in this paper, the transfer function poles 
and zeros of the approximations to the atmospheric 
disturbances can be directly computed for different parameters 
describing atmospheric turbulence. The model is derived to 
duplicate the fractional order form of atmospheric turbulence 
and is more accurate than previous models developed. These 
new models could be used to design controls for aerospace 
vehicle propulsion or flight control systems. 
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Appendix A. — Nomenclature 


a speed of sound, (m/sec) 

a 0 ambient speed of sound, (m/sec) 

c p specific heat capacity at constant pressure, (J/(kg*K)) 

C t equivalent electric circuit capacitance for t-type of 
disturbance, (farads) 

/ frequency, (Hz) 

f n natural frequency, (Hz) 

f c frequency of computed correction factor, (Hz) 

h altitude, (km) 

H vector of intersection points to match estimated and 

fractional order TFs 

k wavenumber, (rad/m) or (cycles/m) 

K proportional gains 

K t von Karman horizontal asymptote for t-type of 
disturbance 

K m adjustment factor to disturbance natural frequency 
/ length, (m) 

L integral length scale, (m) 

L 0 outer scale length, (m) 

m p number of poles in the TF approximation 

m z number of zeros in the TF approximation 

M mach number 

n number of frequency decades desired to estimate the 

fractional order TF 

P pressure, (Pa) 

P 0 freestream static pressure, (Pa) 

P st standard reference pressure, (100,000 Pa) 

q fractional order of equivalent electrical circuit 

R universal gas constant, 287 (N*m)/(kg*K) 

r units conversion exponent of atmospheric disturbance 

R t equivalent electric circuit resistance for t-type of 

disturbance, (ohms) 

S atmospheric turbulence spectral density 

Si- 0 r-v atmospheric velocity turbulence spectral density, 

(m 3 /sec 2 )/(rad) or (m/sec) 3 /(Hz) also converted to 
(m/sec)/(Hz) by taking 1/3 root 

s Laplace operator 

T Potential temperature, (K) 

T 0 freestream static temperature, (K) 

T s static temperature, (K) 

T 0 h temperature at an arbitrary reference height, (K) 
v flow velocity, (m/sec) 

v} velocity fluctuation in a region of size / of the basic 
laminar flow, (m/sec) 


W unity disturbance time domain signal 

W t disturbance time domain signal for t-type of 

disturbance 

x fractional order of atmospheric disturbance spectral 

A.l Greek 

a constant associated with atmospheric turbulence 

spectral density 

a / longitudinal, (unit less), 

a v transverse or vertical, (unit less) 

clt temperature, (°K 2 sec 2 nr 2 ) 

a p pressure, (Pa 2 sec 2 m -2 ) 

y ratio of specific heats, (y = 1.41) 

c eddy dissipation rate, (m 2 /sec 3 ) 

rj the ratio of each decade interval where a pole or a zero 

will be used to estimate the fractional order TF 

X wave length of atmospheric disturbance (m/cycle) 

v kinematic viscosity of a laminar flow, (m 2 /sec) 

p weight density, (kg/m 3 ) 

p pz density of pole-zero pairs per decade for estimated TF 
oo n natural frequency, (rad/sec) 

(Dpt frequency of TF pole i of integer order approx, 
(rad/sec) 

co z * frequency of TF zero i of integer order approx, 
(rad/sec) 

A.2 Subscripts 

a constant associated with circuit capacitance 

correctional factor 

A variable associated with an adjustment 

C constant associated with circuit capacitance 

correctional factor 

H associated with symmetry frequencies for the approx, 
fit variable associated with TF fit or approximation 

/ variable associated with longitudinal atmospheric 

disturbance 

LA variable associated with longitudinal acoustic 
disturbance 

VA variable associated with vertical or transverse acoustic 
disturbance 

P variable associated with pressure disturbance 
t variable associated with type of disturbance 

T variable associated with temperature disturbance 
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TLA 

TVA 

v 

VK 


variable associated with longitudinal acoustic VKA 

disturbance due to temperature 

variable associated with vertical or transverse acoustic o 

disturbance due to temperature p 

variable associated with vertical or transverse z 

atmospheric disturbance 

oo pi 

variable associated with the von Karman spectral 
density °° zz 

00 


variable associated with the von Karman spectral 
density circuit approximation 

variable associated with an output 
variable associated with TF poles 
variable associated with TF zeros 
variable associated with pole frequencies 
variable associated with zero frequencies 
variable associated with freestream conditions 
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Appendix B. — Circuit Analog 

In this appendix the von Karman approximations to the 
different atmospheric disturbances based on the Kolmogorov 
spectrum given by (Soreide and Tank (a) and (b)) are 
approximated by utilizing an electrical circuit analog. The 
purpose of the circuit analog is to serve as basis for deriving 
integer order pole-zero product TF approximations to the 
fractional order atmospheric disturbances. The reason for the 
need to derive approximations to the fractional order equations 
is because of the difficulty of numerically solving fractional 
order differential equations. For instance, unlike integer order, 
fractional order derivatives obey the law of non-locality. This 
property of fractional order derivatives makes the state 
transition matrix a convolution integral, with each integration 
step necessitating the computation of the convolution of the 
state all the way back to time zero where the state had an initial 
value of zero. 


B.l Longitudinal and Transverse Disturbances 

The Kolmogorov atmospheric turbulence spectral model as 
described in the main body, Equation (1), of this paper is 


S t (k) = a t e 2/3 k 5/3 ( B j) 

And the von Karman longitudinal and transverse models 
respectively, reproduced from Equations (5) and ( 6 ) from the 
main body of this report are: 


SiyAk) 


Sv,VK:(k) 


2.7e 2 / 3 Z 5/3 


2.7 s 2/3 7. 5/3 


2 

1 + (l.339(27r)L A:) 2 V6 

l + ^(l.339(27r)7l) 2 
l + (l.339(2Tr)77c) 2 11/6 


(B.2) 


(B.3) 


The shape of the von Karman spectrum approximations of 
Equations (B.2) and (B.3) is that of a parallel network consisting 
of a resistive element corresponding to the horizontal von 
Karman asymptote and capacitive element corresponding to the 
high frequency Kolmogorov asymptote as 

S t ,VKA ( k ) = -^ 7^7 ' , (m/sec)/Hz (B.4) 

{K t y + s[{ k) 

Where r represents the units conversion exponent from 
(m/sec) 3 /Hz to (m/sec)/Hz (i.e., r=l/3 in this case — for the 
longitudinal and transverse acoustic wave disturbances). The 
variable S t (k ) signifies the Kolmogorov spectrum of 
Equation (B.l) or the high-frequency asymptote of the von 
Karman spectrum, and K t> signifies the von Karman low- 
frequency asymptote or TF gain for the type of disturbance of 
Equations (B.2) and (B.3) as 


of the von Karman Model Form 


K, =5.4s 2/3 7, 5/3 
K v = 2.7s 2/3 7 5/3 


(B.5) 

(B. 6 ) 


Equation (B.4), for a flow, such as that of the acoustic velocity 
spectral, is equivalent to a circuit as shown in Figure B.l. The K t 
element is a gain for the type of disturbance and is a function of 
the parameters describing the atmospheric disturbance. The 
parameter W t represents an input flow and W t>0 represents the 
resulting output flow, where the exponent q represents the 
fractional order of these elements. The circuit is an analog of this 
representation, acting as a disturbance source connected to a 
filter, meant to be connected to a sink or a system simulation in 
order to establish an output flow. Applying circuit theory, the TF 
of this circuit can be computed as 


K r 

- ‘ . W t 


(R t C t s) q +\ 


(B.7) 


Where W t in the circuit stands for a unity flow signal 
representing the atmospheric disturbance, containing the 
disturbance frequencies of interest. More on that will be 
discussed later. The output signal, W t>0 , will be the time domain 
signal representing the von Karman spectral approximation. 
The subscript t stands for the type of disturbance as before, the 
exponent q represents the fractional order of the circuit, and the 
factor K t is given by Equation (20) or (21). Where l/(C t s) q and 
R t q are the capacitive and resistive impedances of the circuit, 
with s representing the Laplace operator. 

The values of resistance R h capacitance C h and the fractional 
exponent q remain to be determined. The natural frequency of 
this circuit, co„, can be directly obtained by inspection of the TF 
of Equation (B.7). This frequency will coincide with the 
frequency at which the low-frequency asymptote (Figs. 2 and 
3) intersects the high-frequency asymptote associated with the 
Kolmogorov spectral or the capacitive impedance of the circuit 
as 



R, q 

-vwv 

i/(Qs f 


o 


* W,,o 


Figure B.l . — Equivalent circuit TF of atmospheric 
disturbance. 
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CO; 


(B.14) 


For now K m is used as a placeholder in case any adjustments are 
needed to this natural frequency and its value here is set to 1. 

In Equation (B.l), k stands for the wavenumber and it is 
related to the frequency of the disturbance as 

f = — = / c Ma (B.9) 

X 

For the circuit analogy, the natural frequency of the circuit 
should be set to the same as that for the von Karman spectral of 
Equations (B.2 and B.3). By inspection, the natural frequencies 
of both the von Karman longitudinal and transverse 
disturbances are the same and can be calculated by substituting 
Equation (B.9) for / into Equation (B.2) or (B.3) to obtain the 
relationship 

— = 1.339(2ti)I (B.10) 

co v ' Ma 

Then by substituting the magnitude s for 2nf on the right hand 
side of Equation (B.10) and also setting co„ = 2 nf n on the left 
hand side, we can solve for the natural frequency as 

. Ma 

* ~ 1 .339(271)1 (R11) 

Additionally, the TF of Equation (B.7) needs to represent the 
von Karman approximation of the Kolomogorov model (i.e., 
the Kolmogorov asymptote). However, the Kolomogorov 
model, as a straight line decreasing with frequency in a log 
scale, can be represented in Figure B.l with a capacitive 
impedance. Thus, Equation (B.l) must be this capacitive 
impedance. This capacitive impedance is raised to the r=l/3 
power to account for the units change from (m 3 /(sec 2 cycle)) to 
(m/sec)/Hz. Therefore, taking Equation (B.l) and raising it to 
the power r, replacing k with f/(Ma ) using Equation (B.9), then 
multiplying both numerator and denominator by (27i) (5/3)r , and 
by setting the resulting equation to l/(C t s) q , the following 
expression can be obtained. 


By inspection of Equation (B.12), since |s| = 2nf, then 
q= (5/3 )r or in general q=xr (i.e., with x=5/3 in this case), then 
the capacitance, C h can be computed as 

C t = 7 77 ) (B.13) 

(a t z 2/3 J x (2nMa ) 

Plugging Equations (B.ll) and (B.13) into (B.8), with 
00 = 2 71 /, R t can be solved as 


R t =\.339(2nia t e 2/3 ) Vx L 

Neither the capacitance C t nor the resistance R t depends on the 
exponent r, utilized for units conversion, which could be 
indicative of the fundamental nature of these parameters. 

With this analysis, the circuit elements for the longitudinal 
and transverse von Karman model approximations, pertaining 
to Equation (B.7) have been computed. Figures B.2 and B.3 
show the von Karman longitudinal and transverse spectral of 
Equations (B.2) and (B.3) respectively and their circuit 
approximation of Equation (B.7). These approximations can be 
improved by increasing the magnitude of the circuit 
approximation (K t in Equation (B.7) for the longitudinal wave 
in order to reduce the high frequency asymptote error at the 
lower frequency. 

To do this, an offset dB level is chosen (AdB) so that the 
resultant horizontal asymptote of the circuit approximation is 
not too high above the corresponding horizontal asymptote of 
the von Karman spectral. This can be done by multiplying 
Equation (B.7) by a proportional gain K a , corresponding to this 
AdB level, as K a =10 AdB/2 °. However, doing that will also raise 
the high frequency asymptote of the circuit approximation by 
the same dB level. To counter this affect, the high frequency 
asymptote can be lowered by the same amount. Lowering the 
high frequency asymptote can be done by lowering the 
capacitive impedance, since, the capacitive impedance is 
inversely proportional to the capacitance, see the circuit in 
Figure B.l. Thus, this will also require multiplying the 
capacitance by the same proportional gain, K a . 

For the transverse approximation, illustrated in Figure B.3, it 
is seen that the high frequency asymptotes of the circuit 
approximation runs in parallel with the von Karman Spectral. 
Therefore, the asymptote corresponding to the transverse 



Frequency, Hz 

Figure B.2. — von Karman longitudinal spectral of Equation (B.2) 
and their respective circuit TF approximation of Equation (B.7), 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 
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Figure B.3. — von Karman transverse spectral of Equation (B.3) 
and their respective circuit TF approximation of Equation (B.7) 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


approximation needs to first be raised by doing the reverse of 
what was done to the capacitance in the previous step (i.e., 
multiplying the capacitance by the reciprocal of a certain 
correction factor, K c . Secondly, the magnitude of the transverse 
approximation can be adjusted upwards by a desired level, K a , 
by performing the two steps described for the adjustment of the 
longitudinal approximation. These correction factors can be 
computed as follows. 


K Ua =10 AdB/2 ° 

(desired AdB level may range from 0 to 6 dB) 


(B. 15) 


Kl,c = 1, K VtC = S v y K (f c )- W v , 0 (f c ) , 


(B.16) 


where f c is the frequency where this correction factor is to be 
computed. With these factors computed, Equation (B.13) for 
the capacitance can be modified utilizing Equations (B.15) and 
(B.16). But the approach would be to keep the values of the 
capacitance and resistance generic and instead incorporate these 
adjustment factors directly into Equation (B.7) as 


W tM = 


K x 


t,A 


\l/q A q 

K t,a 


-Wt 


K tc 
V l ’ c J 


R t C t s 


+ 1 


(B.17) 


where 


K t 


t,A 


- K,, a K[ 


(B-18) 


and Kf (for the type of disturbance) comes from Equations 
(B.5) and (B.6) raised to the power r. 



Figure B.4. — Von Karman longitudinal acoustic wave spectral of 
Equation (B.2) and their circuit TF approximations with different 
magnitudes of adjustment (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Figure B.5. — Von Karman transverse acoustic wave spectral of 
Equation (B.3) and their circuit TF approximations with different 
magnitudes of adjustment (s=8.6e-5 m 2 /sec 3 , L=7 62 m). 


Figures B.4 and B.5 show the longitudinal and transverse 
acoustic wave spectral for the von Karman forms of Equations 
(B.2) and (B.3) and the equivalent circuit adjustments based on 
Equation (B.17) for different magnitude adjustments. The K t>a 
values for these adjustments are based on Equation (B. 1 5). The K v>c 
value used is 1 .424. As discussed before, the important frequencies 
to match are the higher frequencies, above 1 Hz, because a typical 
control system design should have no problem attenuating 
disturbances at the lower frequencies. Therefore, the 6 dB or even 
the 3 dB, adjustments would normally provide sufficient 
approximations to the von Karman spectral. 
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B.2 Temperature Disturbance 

Unlike the von Karman model forms of Equations (B.2) and 
(B.3) for the acoustic wave given by Soreide and Tank (1996), 
no comparable models exist for temperature, temperature 
generated acoustic wave, pressure, or density disturbances. 
Thus the approach here is to utilize the Kolmogorov spectral of 
Equation (B.l) for the acoustic wave disturbance and compute 
the differences from that and those of the Kolmogorov Spectral 
of Equation (B.l) for temperature. Then scale the respective 
von Karman models of Equations (B.2) and (B.3) by these 
differences to generate the von Karman models for temperature, 
pressure and density. 

Figure B.6 shows this scaling, with the resulting von Karman 
spectral for the temperature disturbance. This scaling can be 
done by graphically finding or by computing the difference in 
dB magnitude (at any frequency using Equation (B.l)) between 
the Kolmogorov temperature and the Kolmogorov longitudinal 
spectral, and then multiplying the von Karman longitudinal 
spectral by this scale factor (corresponding to the difference in 
dB magnitude) to come up with the von Karman temperature 
spectral. As a result, the von Karman temperature spectral is 
computed as 


$T,VK 


(k) 


= 7.0s 2/3 Z, 5/3 

[l + (l . 339(271 )L A:) 2 
((k 2 *m/sec)/Hz) 


5/6 ’ 


(B.19) 


The horizontal asymptote for this spectral, based on Equation 
(B.19), is 


K t - 14.0s 2/3 Z 5/3 (B.20) 


The circuit approximation for this model is similar to that of 
Equation (B.17), except for some differences in the 
proportional gains and the fractional exponent as 


Wi 


K 


T,A 


T,oA 


-W T 


f K V 7 * 

K T,a 
K T ,c 


RjCt^ 


+ 1 


(B.21) 


Where K T , a is the same as Equation (B.15) for the desired 
adjustment and 


Kt,a = Kt,u^t ■> k t,c = 1 (B.22) 

As before, the fractional exponent is q=xr. However, this time 
for the temperature disturbance, the units conversion exponent 
is r= 1/2 (see Fig. 1 for the temperature units). As such, the 


fractional exponent for temperature is q= 5/6, instead of q= 5/9 
for the acoustic wave disturbances. Figure B.7 shows a plot of 
the circuit approximation for temperature based on Equation 
(B.21), with K T , a = 1.4125 (i.e., raising the circuit approximation 
magnitude by 3 dB) , compared to its scaled von Karman 
spectral of Equation (B.19). 



Figure B.6. — Kolmogorov spectral of the longitudinal wave and 
temperature and scaling of the longitudinal von Karman 
spectral to come up with the von Karman temperature model 
form (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Frequency, Hz 

Figure B.7. — Kolmogorov and von Karman spectral for 
temperature and its von Karman circuit approximation 
(e = 8.6e-5 m 2 /sec 3 , L = 762 m). 
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By perturbing the relation of the speed of sound and 
temperature, a = -yJyRT , and substituting the resulting 

expression for A a into the Mach number equation of M=v/a , the 
relation between the change in temperature and acoustic 
velocity can be obtained as 


Av = 


MyR . rj-i 

— - — AT 
2 a 0 


(B.23) 


Where y is the ratio of specific heats (~ 1.41), R is the universal 
gas constant (-287 -N*m/(kg*K)), and a 0 is the standard speed 
of sound at that altitude (~ 295 m/sec at supersonic cruise 
altitude). Figures B.8 and B.9 show the von Karman spectral for 
the longitudinal and transverse acoustic waves respectively, the 
temperature gust based on Equation (B.19) turned to either 
longitudinal or vertical acoustic wave velocity gust, and the 
combination of the two. As shown in these figures, for the 
combined acoustic wave the acoustic wave velocity due to 
temperature gust dominates in the low frequency range. At 
higher frequencies, both affect the total and the resultant slope 
is more of a combination of the two disturbances. The question 
in encountering atmospheric turbulence is if it’s possible that 
the two effects (i.e., the acoustic wave due to pure wind gusts 
and acoustic wave disturbance due to a temperature fluctuation) 
can combine at the same time to produce worst case conditions. 
As for longitudinal and transverse disturbances combining, in 
terms of the propulsion system, the assumption is that the 
longitudinal disturbance would be the worst case (i.e., the 
longitudinal and the transverse or vertical would not combine 
at the same time). This is especially true for propulsion systems 
located under the wing, since, the vertical acoustic wave loses 
some velocity when it’s turned to a longitudinal acoustic gust 
via the vehicle wing forebody (Ahsun (2004)). 


B.3 Pressure Disturbance 

For an atmospheric pressure disturbance, as discussed in the 
previous section, the Kolmogorov spectral of Equation (B.l) 
will be used to come up with the von Karman model form by 
scaling the respective models of the longitudinal acoustic wave 
disturbance. In this case, however, the constant a p in 
Equation (1) or Equation (B.l) has a dependency on freestream 
(static) atmospheric conditions. Therefore, the scaling factor 
CLp/ai will be explicitly used to scale Equation (B.2) to come up 
with the von Karma spectral, recognizing that the difference in 


magnitude between the Kolmogorov spectral of the acoustic 
longitudinal velocity and pressure is exactly equal to this factor. 
As such, the von Karman model form for pressure can be 
expressed as 


Sp,VK W — 


2.7(« P /a / )^ 2/3 L 5/3 
((pa 2 *m/sec)/Hz) 


2 

[1 + (l. 339(2 7r)Lk) 2 ] 5 16 ? 


(B.24) 


The horizontal asymptote for this spectral, based on 
Equation (B.24), is 

K P =5A(a P /a,)£ 2/3 L 5 ' 3 (B.25) 


Similar to the Temperature spectral, the circuit 
approximation for this model is the same as that of Equation 
(B.17), except for some differences in the proportional gains as 


W) 


PM 


K 

f v \ l/( * 
K P,a 

Kp, c 


P,A 


-Wp 


(B.26) 


RpCps 


+ 1 



Figure B.8. — Von Karman longitudinal acoustic wave spectral, 
temperature disturbance turned to acoustic velocity, and the 
two combined. 
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Figure B.9. — Von Karman transverse acoustic wave spectral, 
temperature disturbance turned to acoustic velocity, and the 
two combined. 



frequency, Hz 

Figure B.10. — Kolmogorov and von Karman spectral for 

pressure and its von Karman circuit approximation, (s = 8.6e-5 
m 2 /sec 3 , L = 762 m). 

Where Kp >a is the same as Equation (B.15) for the desired 
adjustment and 


Kp,A ~ Kp, a K r p , K P c = 1 (B.27) 

As with the temperature in the previous section, the fractional 
exponent is q=xr , with r=l/2. As such, the fractional exponent 
for pressure is the same as that for temperature (i.e., q= 5/6,). 
With the value of Kp >a = 1.4125. Figure B.10 shows pressure 
comparisons of the Kolmogorov spectral of Equation (B.l), the 
respective scaled von Karman spectral of Equation (B.24) (both 
converted from Pa 2 to Pa), and the respective circuit 
approximation of Equation (B.26). 



Frequency, Hz 

Figure B.l 1 . — Kolmogorov and von Karman spectral for 
density and its von Karman circuit approximation (e = 8.6e-5 
m 2 /sec 3 , L = 762 m). 


B.4 Density Disturbance 

The density disturbance won’t be needed when acoustic 
velocity, temperature and pressure disturbances are included in 
a simulation. However for completeness, density disturbance 
will be discussed here in brief. The magnitude of the density 
disturbance can be computed from the equation of state, 
P=pRT, by perturbing this equation, which can be solved as 


Ap = 


p 0 + ap 

R(T 0 + AT) 


Po_ 

RT 0 


(B.28) 


where P 0 and T 0 are the standard atmospheric pressure and 
temperature in Pa and Kelvin respectively, which would be 
around 5500 Pa and 216 K at supersonic cruise. Figure B.ll 
shows a plot of the Kolmogorov, the Von Karman, and the 
circuit approximation spectral for density, all based on Equation 
(B.29) from values of temperature and pressure obtained in the 
previous sections. It can be seen from this figure that the 
Kolmogorov spectral obtained from Equation (B.28), actually 
looks like a finite von Karman spectral at low frequencies. Also, 
its low frequency asymptote would be approximately -20 dB, 
which equates to a density of 0.1 kg/m 3 . From standard 
atmospheric tables, the atmospheric density at about 62,000 ft 
(i.e., at approximate cruise altitude) is also approximately 0.1 
kg/m 3 . Even though, at the limit as frequency approaches zero, 
the Kolmogorov spectral approximates the mean atmospheric 
density at cruise, the appropriate density disturbance would be 
that shown for the von Karman or the circuit approximation. 
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Appendix C. — Fractional Order TF Fits of Atmospheric Disturbances 


Normally, with the derivation of circuit analogs and their 
respective TFs to describe atmospheric disturbances (as they 
were presented in the Appendix B), time domain simulations 
would be straightforward. However, these circuits and their TFs 
are fractional, which complicates performing time-domain 
simulations. The complication comes from a certain property of 
fractional order differential equations, which deals with non- 
locality. There are several works that deal with time domain 
solutions of fractional order differential equations Schmidt 
(2006) and Lorenzo (2008), but even to date, solving these 
types of differential equations is still challenging. Unlike 
integer order differential equations, the solution space of 
fractional order derivatives is not local and it depends on the 
whole prior history of the state, resulting in a state transition 
matrix that is in the form of a convolution integral. Therefore, 
solving the convolution at each time instant, taking into account 
all the prior history of the state, can be both complicated and 
taxing on the computing resources. 

In this section the attempt will be to formulate an integer 
order approximation of the fractional order TF. Given the 
parameters describing an atmospheric disturbance, like eddy 
dissipation rate and integral length scale, the formulations will 
explicitly solve for the integer order TF poles and zeros that 
approximate the fractional order disturbance. This approach 
avoids the laborious process of hand fitting such an 
approximation, for every type of disturbance, and for every time 
an atmospheric disturbance parameter changes. After forming 
this approximation to the fractional order TF, and by selecting 
the frequencies of interest for the disturbance as an input to this 
TF approximation, the time domain atmospheric disturbance 
can be formed. 

If the desire is to estimate the fractional order TF at each 
frequency decade by a single first order pole-zero pair, then the 
estimation process for each frequency decade starts with a pole 
followed by a zero (i.e., co /; /< oo z *), with one last pole placed after 
the last decade in order to keep attenuating the TF magnitude as 


W t , 0 =K 


f „ \ 

TT s ! ®zi + 1 

. i S/CDjyi + 1 

\l=l / P l 


/®pn+\ 


+ 1 


w u I = 1,2 /t (C.l) 



Figure C.l. — Pictorial diagram of fractional order TF fit. 


would be decreasing at a rate g*20 dB/decade, where q is the 
fractional order. Therefore, midway through the first decade the 
initial gain is expected to decrease by a factor of (%)*</ *20 dB 
as 


K 


s/topi +1 


—q{ 20)/20 

: 10 2 ^ K 


(C.2) 


Divide by 20 is because of the 201ogio(x) base scale. Solving for 
the frequency we get 


I _ s 

°M,=l/2decade JQ(1/2)^_| 


(C.3) 


Also, midway in the first decade, the fractional TF, based on its 
circuit approximation in Appendix B, would have the same gain 
as that of Equation (C.2) as 


K 

(R t C t sY+\ 


.s-l / 2 decade 


-1 

= 102^ 


(C.4) 


As was discussed in Appendix B, the natural frequency of the 
fractional order TF in Equation (C.4) is 


Where n is the number of pole-zero pairs of the estimated TF 
and K is the DC gain. So if the desire is to approximate the TF 
for three consecutive decades, with one pair of pole/zeros for 
each decade, starting at some desired frequency, then n would 
be three in this case. A sketch of this TF approximation as a 
staircase symmetrically located on top of the fractional order 
TF is shown in Figure C.l, with more description to follow 
later. 

Thus, if one first order pole-zero pair is used to estimate each 
decade of the fractional order TF, the first pole should be placed 
at a certain frequency such that the pole gain dropping at a rate 
of -20 dB/decade, intersects the gain of the fractional order TF 
midway through, in the first decade. A fractional order TF gain 


Knowing co„, Equation (C.4) can be solved for s to find out at 
what frequency to place the first pole in order for the gains of 
these two TFs to intersect midway through the first decade as 

M =co„(lO< 1/2 )<?-lT /<? (C.6) 

5=1/ 2 decade 

Substituting Equation (C.6) into Equation (C.3) for s, the 
frequency of the first pole can be solved as 
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(C. 12) 


<V 


= CO„(l 0(1/2), _1 jf 


(C.7) 


^ /( ° zl+1 ^ | |d d = 10 r{H*)q K 
s/(to p \+l 1 


So far in this development the assumption was that each decade 
of the fractional order TF will be approximated with one pole- 
zero pair as shown in Equation (C.l). Based on that, the 
estimated product of first order poles and zeros of Equation 
(C.l) will match the fractional order TF at half decade intervals. 
For that, the only distinction in Equation (C.7) is the 1/2 
exponent in the right side of the equation. This successive ratio 
of desired decay matching, expressed here as //, can be 
generalized as 


<o p1 =o)„(lO (w ^-l)T (C.8) 

where H p \ is the decay ratio for the first pole, where the two TFs 
would match or the magnitudes will intercept. 

Following the placement of the first pole, the need is to find 
out at what frequency to place the first zero of this pair. At the 
end of the first decade, the gain of the fractional TF would have 
dropped by 10 ~ q K as 

7 V/ s at first decade —10 ^ K (C.9) 

{ R tC t s) q + 1 

By substituting 1/c o/, using Equation (C.5), for R t C t , the 
frequency s where the gain has dropped by the above magnitude 
can be calculated similarly as Equation (C.6) 

l.vl =co„(l0?-l) 1/ ‘ ? (C.10) 

at 1 decade 


Or in general terms, as a function of the second portion of the 
approximation, H z \, desired to be matched (i.e., the first portion 
at Vi decade was matched with a pole, the second portion at 1 
decade is matched with a zero (H z i=l), and so on), this 
frequency can be expressed as 

®Hzi = (l 0 (//rl )q - 1 )' ' q (C.ll) 


where 


From this equation, the frequency of the first zero can be 
calculated as 


®*i = 


®Hz\ 

10 ^ zl ^(co// z i/co^i +l)-l 


(C. 13) 


Thus, choosing the ratio of decades at which it’s desirable to 
match the TF estimate, Equations (C.2), (C.8), (C.ll), and 
(C.13) can be used to compute the frequency of the first zero. 

Following the same procedure to calculate the frequency of 
the second pole, by using Equation (C.9), but with the gain 
10-( 3/ 2)^, (i.e., H p 2=1.5 decades), 

^H P 2=^M Hpl)q (C.l 4) 


Then using the TF with two poles and one zero equated to this 
gain, similar to Equation (C.12), the frequency of the second 
pole can be calculated as 


® p2 ~ 


^Hpli^Hpl /topi + 1) 

I ()( h p2 h (m Hp 2 /co-i + 1 ) —1 


(C-15) 


The same procedure can be followed to calculate the rest of the 
frequencies of the poles and zeros of the estimating TF. In order 
to generalize these formulations (for approximating the 
fractional order TF), let’s define p pz as the desired density of 
pole-zero pairs per frequency decade. In addition, let’s also 
define the number of equal decade subintervals (in log 
frequency scale) as r|, where a pole or a zero will be used to 
closely match this TF as 


ti = -E (C.l 6) 

2 P pz 

Then, pole and zero vectors of decade intervals can be formed, 
where these TF will be closely matched as 


(toHzl n^// zl (decade) 

The purpose of the cdh frequencies is to ensure symmetry (i.e., 
to symmetrically locate the TF approximation on top of the 
fractional order TF. Thus, in order to accomplish symmetry, the 
desire is to also have the estimating TF with the first pole-zero 
pair to equate to the same gain at this frequency, oohzi, (i.e., the 
same gain as the right hand side of Equation (C.9), see also 
Figure C.l, as 


H p = r| [2(1) - 1 2(2) - 1 . . .2 (m p ) - 1] (C. 17) 

H Z = \\\2(X) 2(2). . .2(m z )\ (C.18) 

For n number of decades over which to estimate the fractional 
order TF, m p and m z can be computed as 

m p = (n— l)/r| (C.l 9) 

m z = m p - 1 (C.20) 
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For instance, for a density of one pair of pole-zeros per decade 
(i.e., p pz =1) and estimating the TF for three decades (i.e., n= 3), 
these vectors would have the values H pi = [1/2 3/2 5/2 7/2], H zi 
= [12 3], 

Inspecting Equation (C.17), it can be seen that the first 
element of the vector, corresponding to the calculation of the 
first pole frequency, will always be equal to the value of rj. 
Therefore, for convenience, Equation (C.8) can also be 
expressed in terms of r\. 

® pl =<s> n (lOW-lYr (C.21) 

Also, inspection of the equations that describe the frequencies 
of the zeros and poles (i.e., Eqs. (C.ll), (C.13), (C.14) and 
(C.15) respectively, it can be seen that a generalized 
relationship for these equations can be written as follows 

oo^=co„(lO^'V-l)f /? (C.22) 


®zi 

Kg zi^Hzi (© Hzi / m z/-1 + 1 \®Hzi / ©zi-2 + 1 / ©z/1 + Q 

10 ( H zi)q [(cD^/cD^- + iJtoHzi/topi-l + 1 + l)]“ 1 

l-i (C.27) 

K(£>zi © Hzi | | Hzi /© zi-j 0 

= — for i = l,2,...m z 

io- 2 ^Yl^ Hzi /(D Pi +i)-i 

7=1 

These formulations complete the computation of the pole-zero 
frequencies for the TF approximation of fractional order 
atmospheric disturbances. In terms of the indexes defined here, 
Equation (C.l) can be rewritten as 

m z 

fl( S /®ZJ +1 ) 

W t ,o = K tJlt -J- W t (C.28) 

If(V^ + l) 

1 


A more generalized expression of Equation (C.22) can be 
arrived by substituting into it H p from Equation (C.17) in terms 
of the pole number, i and r\ as 

®H P i=®n( 10 1 i (2 *'- 1) ‘?-l) 1/ ‘? i = 2,...,m p (C.23) 

Then a generalized expression for the calculation of the pole 
placement frequencies can be formulated as 


As described in Appendix B, the von Karman spectral, 
Equations (B.2) and (B.3), can be approximated with a 
fractional order circuit with a natural frequency that can be 
determined using Equation (B.8). Equation (B.8) is 
conveniently reproduced below in Equation (C.29). 


= 

_ ^api^Hpi ip^Hpi pi- 1 ^fp^Hpi j © pi-2 /©/?il + l) 

10 ( H pi) q ]^co Hpi / CO zi _i + l)(a > Hpi / ©zi- 2 + l)....(©// p / / ©zl + 1)]— 1 
i-l 

7=1 . 0 

= — x . i = 

io^n(^-Kw + i)-i 

7=1 


(C.24) 


The proportionality factors Koyi have been inserted in case any 
final adjustments are needed to these frequencies to improve the 
TF fit. 

Similarly, for the frequencies for placing the zeros 

® 7 /z, =co„(l0 ,w --'> / -l) l/ ' ? i = \,2,...,m z (C.25) 

Or by substituting in Equation (C.21) the H z vector from 
Equation (C.l 8) in terms of the zero number, /, and r| as 

©/&,•=©„(] l0 2 ^-lj lq i = l,2,...,m z (C.26) 

Then a generalized expression for the calculation of the zero 
placement frequencies can be formulated as 


with an equivalent capacitance and resistance for the type of 
disturbance t , derived as Equation (B.13) and (B.14) 


(a t E 2/3 J /x (2nMa) 

(C.30) 

= 1.339(2ji )(a t s 2 / 3 f x L 

(C.31) 


Based on these derivations, TF approximations of fractional 
order atmospheric turbulences can be developed (for the von 
Karman spectral, for any fractional order 0<<?< 1) by using 
Equation (C.28), together with supporting Equations (C.l 6), 
(C.l 9) to (C.21), (C.25) to (C.27), (C.23) to (C.24) and (C.26) 
to (C.27), (C.29) to (C.31). 

In a pictorial sense, by referring back to Figure C.l, the first 
order pole-zero pair TF approximation of a fractional order 
atmospheric turbulence can be viewed as a staircase 
symmetrically located on top of the fractional order TF. The 
utility of (Hi i P i and oo Hzi are to maintain this symmetry. This is 
accomplished in the derivations above by formulating the 
symmetry frequencies, oo/fs, such that the length segments 
between intercepts (like d 1 and d2) or alternatively the 
horizontal distances on either side of the intercepts are equal. 

In the subsections that follow, more detail will be discussed 
about developing the approximations for each type of 
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disturbance. In addition, some final adjustments that need to be 
made to the pole-zero frequencies and the proportionality 
constants K will be presented. Without these final adjustments, 
the fits will resemble the example fit shown in Figure C.2, with 
the fractional order fit winding its way around the fractional 
order TF 


C.l Longitudinal and Transverse Acoustic 
Wave Disturbances — TF Fits 


By applying the formulations derived in Section C.l, TF fits 
for the longitudinal and transverse atmospheric disturbances 
can be constructed based on Equation (C.2 8). Figure C.3 shows 
a plot of the longitudinal von Karman disturbance of Equation 
(B.2), it’s circuit TF approximation reproduced from Figure 
B.2, and the TF fit computed using the equations derived 
previously in this Appendix for q=xr= 5/9 (i.e., x=5/3, r= 1/3), 
r] = l/2, n= 3 (decades), and with atmospheric turbulence 
parameters s=8.6e-5 (m 2 /sec 3 ) andT=762 m. 

The maximum error in this TF fit for the circuit 
approximation, using 4 poles and 3 zeros, over a span of VA 
decades, amounts to approximately 1.5 dB. The error for the 
actual von Karman spectral is larger; about 4 dB at 0.1 Hz and 
decreases as the frequency increases. Assuming, that the control 
system design can sufficiently attenuate frequencies at this low 
frequency range (as it should) then this TF fit approximation 
can be acceptable for the actual von Karman spectral. For the 
TF fit shown in Figure C.3, which is an approximation to the 
longitudinal von Karman spectral of Equation (B.2) the pole 
and zero frequencies in Equation (C.2 8) are 


(Dpi = [0.6 12.54 85.71 2062.0], 

r i (rad/sec) 

co zi ~ [3.82 22.92 312.38] 


(C.32) 


with Ki t f lt in Equation (C.28) given by the numerator in 
Equation (5) raised to 1/3 power to account for the units 
conversion as 


Ki, fit = 


(5.46 2/ V 5/3 




(C.33) 



- 2-1012 
10 10 10 10 10 

Frequency, Hz 


Figure C.2. — An example Fit without final adjustments to 
proportionality factors and frequencies. 



Frequency, Hz 

Figure C.3. — Longitudinal von Karman spectral, circuit 
approximation, and TF fit (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


The preceding development to derive the equations to 
approximate a fractional order TF was based on a single 
fractional order TF, and a circuit approximation to the 
atmospheric disturbance was used for this development. 
However, as indicated by Equation (B.2) and (B.3), the von 
Karman spectral is not that of a single order fractional TF, 
especially in the lower frequency spectral. Therefore, for better 
accuracy, the derived equations can be adjusted to better fit the 
von Karman spectral. Examining Figure C.3, it can be seen that 
the TF fit can better approximate the von Karman spectral if the 


natural frequency oo„, Equation (C.5), is increased, about 
proportionally to the difference in magnitudes between the 
fractional order TF and the approximation TF, at that frequency. 
However, doing that will also increase the pole-zero 
frequencies of the TF fit based on the equations derived. 
Therefore, some equivalent adjustment to some of these 
frequencies will need to be made. Figures C.4 and C.5 show the 
TF fit for the longitudinal and transverse von Karman 
disturbances respectively, with the following adjustments to 
Equations (C.24), (C.29) and (C.27) 
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K, 


7/v, co 


\K(jM 5 Kfopi , J 

[2.4; 1 1 1/2.4 1/1.5; 1 1 l] 


(C.34) 


With Ki/ V , a symbolizing the proportionality constants for the 
frequency adjustments for both the longitudinal and transverse 
disturbances, which gives the frequency values of ( for s=8.6e- 
5 m 2 /sec 3 andT=762 m) 


(a pi = [l .46 30.10 85.71 1593.1], 
co z! =[9.18 55.02 335.48] 


(C.35) 


The proportionality factor, Ki t f lt , in Equation (C.28) is given 
again by Equation (C.33). Similarly, K v> f lt is developed using 
Equation (B.6) as follows, but also adjusting by 3 dB which 



Figure C.4. — Longitudinal von Karman spectral, adjusted TF fit 
(e = 8.6e-5 m 2 /sec 3 , L = 762 m). 


will appropriately raise its magnitude in order to match the high 
frequency asymptote as 

K v ,fit = l-4(2. 7s 2/3 I 5/3 )l (C.36) 

As can be seen in Figure C.5, the TF fit of the transverse 
disturbance is not very accurate at the low frequency, which 
would normally be fine for control system design purposes as 
discussed before. But it may not be acceptable, for let’s say, wing 
loading. Thus, further adjustments could be made if desirable. 
For instance, instead of multiplying K v> f lt with the 3 dB gain of 
1.4, the natural frequency oo„ of Equation (C.29) could be 
multiplied by this factor to raise its frequency. Then for final 
adjustments, the natural frequency can be also multiplied by the 
ratio of the difference of frequencies between the fit and the 
actual von Karman spectral obtained at some fixed amplitude, at 
some high frequency. This can be done either graphically or 
analytically. Performing these manipulations, a better fit of the 
transverse disturbance is shown in Figure C.6, with a final value 
of K m = 4.27. This will result in the following adjustments for 
the transverse disturbance 

K v ,a — \K om , K 0}pi , K^i ] ^ ^ 

= [4.27; 1 1 1/2.4 1/1.5; 1 1 l] 

which gives the frequency values of ( for s=8.6e-5 m 2 /sec 3 and 
L=1 62 m) 


(o V9 pi = [2.60 53.56 152.55 2835.3], 
(o VZ i = [l6.33 97.92 597.07] 

with a K v> fit expression as 

*v.fi t =( 2.7s 2/3 Z 5/3 )i 


(C.38) 


(C.39) 




Figure C.5. — Transverse von Karman spectral, adjusted TF fit Figure C.6. — Transverse von Karman spectral, final adjusted 
(e = 8.6e-5 m 2 /sec 3 , L = 762 m). TF fit (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 
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C.2 Temperature Disturbance — TF Fit 

As discussed before, the units conversion factor for 
temperature is r= 1/2, which makes the fractional exponent 
value g=5/6 (*=5/3). Substituting this value for the fractional 
exponent, with the same values of rj, n , and with the same 
atmospheric parameter values as those used for the longitudinal 
and transverse acoustic disturbances, and again using the 
formulations previously developed in this Appendix, the poles 
and zeros for the TF fit can be directly computed as before. The 
resulting TF fit with K t> f lt in Equation (C.28) given by Equation 
(B.20) raised to the Vi power as 

K t =(l4.0e 2/3 Z 5/3 )2 (C.40) 

This fit (not shown here) will end-up with the same low and 
high frequency asymptotes as with the von Karman form, but at 
somewhat lower frequency. Similarly to the longitudinal 
disturbance, this can be corrected by multiplying K® n in 
Equation (C.29) by 3/2, which turns out to be the ratio of the 
frequencies of these two spectra at some high frequency, for 
select fixed amplitude. Finally, by some minor adjustments to 
the high frequency poles, these spectra can be matched fairly 
close, as shown in Figure C.7. This will result in the following 
adjustments for the temperature disturbance 

K-T , cd — \Kq on 5 ^cq pi ? ozi ] (C 4 1 ) 

= [l.5; 1 1 1/El 1/1.2; 111] 

which gives the frequency values of (for s=8.6e-5 m 2 /sec 3 and 
L = 762 m) 


©^•=[1.10 25.11 109.77 816.35], 
CD t zi = [33.04 45.64 602.36] 


(C.42) 


with a Kt, fit, using Equation (B.20) as 

Kt, fitctemp) = (l4.0£ 2/3 L 5/3 y /2 (C.43) 


If the desire is to simulate the acoustic disturbance generated 
by a temperature gust, this conversion can be done by utilizing 
the perturbation relation of the speed of sound with temperature 
described in Appendix B (Eq. (B.23)) as 


MyR „ 
Av = — —AT 

2a 0 


(C.44) 


Together with Equation (C.36) to adjust Kt, fit as 

K T , fit(acoustc) = ^ Vl4.0 £ 2/3z5/3 ( C .45) 

AG,q 



Figure C.7. — Temperature von Karman spectral and its TF fit 
(s=8.6e-5 m 2 /sec 3 , L=762 m). 



Frequency, Hz 

Figure C.8. — Von Karman acoustic wave velocity spectral due to 
temperature gust and its TF fit (s = 8.6e-5 m 2 /sec 3 , L = 762 m). 


The acoustic velocity disturbance due to a temperature gust, at 
approximate supersonic cruise altitude is shown in Figure C.8. 

Based on all the derivations carried out so far, the TF fits 
should work for different values of atmospheric parameters 
(like those of Eq. (C.45)). Figures C.9 and C.10 as compared to 
Figures 3 and 4 for their shape, show the TF fits with different 
values of eddy dissipation rates and integral scale lengths, by 
employing the same TF fit equations and without changing any 
other parameters, like the Kt,® values, Equation (C.41). 
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Figure C.9. — Von Karman acoustic wave velocity spectral due 
to temperature gust and its TF fits for different integral scale 
lengths (e = 8.6e-5 m 2 /sec 3 ). 



Frequency, Hz 

Figure C.1 1 . — Pressure von Karman spectral and its TF fit 
(s = 8.6e-5 m 2 /sec 3 , L = 762 m). 



Frequency, Hz 

Figure C.10. — Von Karman acoustic wave velocity spectral due 
to temperature gust and its TF fits for different eddy 
dissipation rates (L=7 62 m). 


(C.27), (C.23) to (C.24) and (C.26) to (C.27), (C.29) to (C.31), 
the poles and zeros for the TF fit for pressure can be computed 
as before. The resulting TF fit, with K t> f lt in Equation (C.28) 
given by Equation (B.25) raised to the V 2 power, is shown in 
Figure C. 1 1 . It turns out, that adjustments to Kp jG) are exactly the 
same as those for Kp,®, Equation (C.41). For completeness, 
these frequency adjustments and the Kp t f lt equation for pressure 
disturbances are listed below. 

Kp , co — [K om , K^pi , ^ ^ 

= [l .5; 1 1 1/1.1 1/1.2; 1 1 l] 

which gives the frequency values of (for c=8.6e-5 m 2 /sec 3 and 
L=162 m) 


(Op p i = [l.lO 25.11 109.77 816.35], 
CO pzi = [33.04 45.64 602.36] 


(C.47) 


with a Kpf lt , using Equation (B.25) as 

Kp i& = Vll.6e 2/3 I 5/3 (C.48) 


C.3 Pressure Disturbance — TF Fit 

The units conversion factor for pressure is the same as that 
for temperature (i.e., r=l/2), which also makes the fractional 
exponent value q=5/6. Substituting this value for the fractional 
exponent, with the same values of r|, n , and //as those used for 
the longitudinal, transverse and temperature disturbances, and 
again using Equations (C.16), (C. 19) to (C.21), (C.25) to 


C.4 Density Disturbance — TF Fit 

As discussed in Section B.4, once the temperature and 
pressure TF fits have been calculated based on the formulations 
in the previous two sections, the density disturbance can be 
directly obtained by utilizing Equation (B.28). A plot of this 
density spectral is shown in Figure C.12. 
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Frequency, Hz 

Figure C.12. — Density von Karman spectral and its TF fit 
(e = 8.6e-5 m 2 /sec 3 , L = 762 m). 
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